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Introduction 

Research in aeroelasticity with computational methods, including the nu- 
merical prediction of flutter boundaries, requires a sophisticated algorithm 
and computer code to simulate unsteady flow phenomena in three dimen- 
sions. The purpose of this research is to develop an improved algorithm and 
computer code for solving the unsteady three-dimensional Euler and Navier- 
Stokes equations. These equations accurately describe the flow phenomena 
for aeroelastic applications. 

An algorithm that uses central differencing has being used in the com- 
puter code, ENSAERO, for aeroelasticity at the NASA Ames Research Cen- 
ter. However, modern upwind algorithms can produce more accurate solu- 
tions for the Euler and Navier-Stokes equations. Among upwind algorithms, 
a streamwise upwind algorithm has recently been developed by Obayashi 
and Goorjian (see Appendix A). Most multidimensional upwind algorithms 
are first constructed in one dimension and then extended to multiple di- 
mensions by applying the one-dimensional procedure in each coordinate di- 
rection. On the other hand, the present method follows the flow physics 
more closely than the coordinate upwind methods. The steady-state com- 
putations confirmed the higher resolution of the present algorithm over the 
centred-difference method as well as over other upwind methods. 

The objective of this research is to update Ames’s aeroelasticity code, 
ENSAERO, by using the improved streamwise upwind algorithm. This re- 
port summarizes briefly the work performed during the period, April 1, 1989 
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through March 30, 1991. Additional details on the various aspects of the 
study are given in Appendices. 


Research Efforts 

The following specific objectives have been performed: 

• The algorithm was extended to handle the moving grid system. The 
finite-volume concept is essential to extend the algorithm. The result- 
ing algorithm is conservative for any motion of the coordinate system. 

• Two extensions to an implicit method were considered. The implicit 
extension that makes the algorithm computationally efficient is imple- 
mented into ENSAERO. 

• The new flow solver has been validated through the solution of test 
problems. The test cases include three-dimensional problems with fixed 
and moving grids. 

• Finally, the new flow solver has been implemented into NASA Ames’s 
aeroelasticity code, ENSAERO. 


Results 

In this work, only the first-order time-accurate methods are considered be- 
cause of computational efficiency. However, time accuracy is an essential 
requirement for aeroelastic computations. Numerical schemes used for flow 
calculations in aeroelasticity must guarantee the correct calculation of am- 
plitude and phase of unsteady pressures. In order to verify the time accuracy 
of the present code, unsteady flows over various wings undergoing prescribed 
oscillatory and ramp motions have been computed. 

F-5 Wing 

The first test case shown here is an unsteady viscous flow over the F-5 wing 
which has an aspect ratio of 2.98, a taper ratio of 0.31 and a leading edge 
sweep angle of 31.92°. The test case is chosen at Af a 0 = 0.9, Tie — 9 x 10 8 
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and the angle of attack of 0°. Since the airfoil section is supercritical, the 
flow is shock-free on the most of the wing surface at this flow condition. In 
the unsteady case, the wing is pitching with a pitch amplitude of 0.11° at 
a reduced frequency of 0.55 as illustrated in Fig. 1. The experimental data 
show that the shock wave appears on the wing surface due to the pitching 
motion. Figure 2 shows the comparison of the computed instantaneous pres- 
sure coefficient contours using the upwind and central-difference methods 
at a = 0° during pitch down motion. The upwind result shows the shock 
wave moving forward, which is captured within the two grid points here. 
In contrast, the central-difference result does not show any concentration of 
contours throughout the cycle. The result confirms that the upwind method 
predicts the shock motion more accurately than the central-difference method 
does. 


Clipped Delta Wing 

The second test considers the motion of the leading-edge vortex as well as 
the motion of the shock wave. Figure 3 illustrate the planform of a clipped 
delta wing and a typical flow field to be computed. Since the leading edge 
is sharp and swept, a leading-edge vortex is formed at moderate angles of 
attack. A shock wave is also formed at transonic speeds. They interact with 
each other at certain flow conditions. 

Figure 4 shows the comparisons of unsteady pressures between the com- 
puted and measured data at M 0 0 = 0.9, Re = 18 x 10 6 , the mean angle of 
attack of 4°, the pitch amplitude of 0.5° and the reduced frequency of 0.6. 
Two peaks in the C v magnitude and the corresponding jumps in the phase 
angle are observed. The first peak near the leading edge is due to the mo- 
tion of the leading-edge vortex and the second peak is due to the motion 
of the shock wave. The results indicate that the computation captures the 
main structure of the flow field with the modified Baldwin-Lomax turbulence 
model. 

To check the aeroelastic option of the code, the flexibility of the wing is 
next considered. Figure 5 shows the first four mode shapes of the clipped 
delta wing. Figure 6 shows the pressure responses on the wing upper surface 
in the 10° ramp motion from 0 angle of attack at M 0 0 = 0.9. At first, the 
leading-edge vortex and shock wave develop. At the inboard section, they 
remain without an interaction. At the outer sections, they have a strong 
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interaction, which leads to a vortex breakdown. Figure 7 shows stream line 
patterns before and after the breakdown. At 1500 time steps ( t = 0.1 sec), 
the angle of attack reaches 10° and the wing is deformed in the nearly maxi- 
mum displacement. Then, the breakdown proceeds from the wing tip to the 
middle of the wing span. At 2400 time steps (t = 0.16 sec), the wing is less 
deformed because of the reduction of the local lift due to the breakdown. 
The results successfully demonstrate the capability of the upwind version of 

ENSAERO. 

For the aeroelastic case, the upwind computation requires 19 fisec per grid 
point per time step at a speed of 160MFLOPS on a CRAY-YMP computer 
using a single processor, while the central-difference computation requires 17 
fisec at a speed of 150 MFLOPS. 

Concluding Remarks 

In the present study, a streamwise upwind algorithm has been extended and 
validated for solving the unsteady three-dimensional Navier-Stokes equa- 
tions. The resulting algorithm has been implemented into Ames’s aeroe- 
lasticity code, ENSAERO. The upwind version leads to higher accuracy in 
both steady and unsteady computations than the previously used central- 
difference method does, while the increase in the computational time is small. 

The future research plan is to continue the algorithm enhancements of 
ENSAERO and to apply the code to complicated geometries, such as a wing- 
body combination and a wing with a control surface. ENSAERO inherits the 
zonal grid capability developed for the Transonic Navier-Stokes code to han- 
dle complicated geometries, such as wing-body and complete aircraft config- 
urations. The key development of the algorithm is integration of the upwind 
and zonal techniques. This includes an interface treatment between moving 
zones for a control surface. The other key development is an improvement 
in turbulence modelling. So far the Baldwin-Lomax model has been used. 
Applicability of more accurate models to unsteady transonic flows is to be 
investigated. 
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Fig. 1 Unsteady modal motion of F-5 wing. 





Circular-arc airfoil 
t/c = 0.06 

L.E. sweep angle = 50.4 R 
Area = 1635.88 in 2 
Span = 45.08 in. 

Root chord = 63.55 in. 
Tip chord = 9.03 in. 
Taper ratio = 0.1421 



Fig. 3 Planform geometry of clipped delta wing and typical flow structure. 



Clipped Delta Wing 

M = 0.90, Re = 18 x 10 6 , a =4° 

°° m 
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□ Experiment, 
NASA TP-2594 



Fig. 4 Comparison of computed unsteady pressures using the standard and 
modified Baldwin-Lomax turbulence models with experiment, = 0.90 
a m = 3.97, d = 0.46, Re c = 17.6 x 10 6 , k= 0.5919. 




Mode 1 : Frequency = 6 Hz 



(b) Mode 2: Frequency = 8 Hz 



( c ) Mode 3: Frequency = 18 Hz 
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Mode 4: Frequency = 28 Hz 


Fig. 5 First four mode shapes and frequencies of clipped delta wing 
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as follows: 
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Numerical Algorithm 


The vector of conserved quantities Q and the inviscid flux 
vector F are 
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where h is the total enthalpy and the contravariant velocity 
component, (7, is defined as U = rf z u 4 tj v v 4 tj x w. For the 
£ and £ directions, E and G can be defined similarly. The 
viscous flux vector G v is given by 


G v = 
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FVS and FDS Algorithm 

The space-discretized form of Eq. (3) can be written as 
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where a second-order central-difference evaluation is applied 
to the viscous term. Now the streamwise upwind algorithm 
is described for the inviscid fluxes at cell interfaces. 

For subsonic flows, Goorjian’s FVS algorithm 4 can be 
written as 




dQ } (6) 


with 

= (l + (j + Ct 

m 2 = 4 

where 72e is the Reynolds number, Pr is the Prandtl num- 
ber, c is the speed of sound, and J is the transformation 
Jacobian. Pressure is related to the conservative flow vari- 
ables, (j, through the equation of state for a perfect gas: 

p = (7 - 1 ){ e - -(u 2 4 v 2 4 w 2 ) } ( 2 ) 

where p is the fluid density and e is total energy per unit of 
volume of the fluid. 

Assuming conical symmetry, Eq. (1) can be rewritten as 
d r Q 4 d v F 4 d<G = -2 E 4 (3) 

where £ is the conical direction and (* 7,0 are the generalized 
coordinates on the conical surface. The conical transforma- 
tion can be described as 9 £ = x A, 77 = tj(Y, Z),( — C(F, Z) 

where A = (1 4 Y 2 4 Z 2 )’, Y = y/z, Z - z/x . The metrics 
are expressed as 



Tj z = -A (Ytjy + Z*iz), Vy = = *nz 

C* = -A(FCv 4 Z(z) t <v = Ky, <x = A<z 


where Q\ and Q r are left and right states, respectively, and 
the metric terms, tj z ^tj v and 7j ZJ are normalized by |Viy| as 
k z = and so on. The contravariant velocity is also 

normalized by IV 77 I, and used as U = k z u 4 k v u 4 k x w in the 
following. For the first-order-accurate computations, l = j 
and r = j 4 1. For higher-order extensions, the MUSCL 
approach 3,8 is used. 

The integral part was derived as follows for U > 0: 

I" 1 ’ ^dQ = {F+A-F} r -{F + A'F}, (7) 

Jqi °q 

where * indicates local sonic values 4 and A m F = F* — F . 
Using the local isentropic relations, A *F can be simpli- 
fied in the one-dimensional case: A *F — A*(pu)e., where 
e # = (\iU t h) T is a sum of two eigenvectors of A*F/A*Q. 
Note A*(pu) > 0 for M < 1 and A m (pu) = 0 for M = 1 . 
This formula is equivalent to the FDS formula, if the flow 
is isentropic. If not, it is an FVS algorithm. Even in the 
latter case, the use of isentropic relations does not restrict 
the flow fields because the isentropic relation is used only 
locally at each grid point (or cell) to compute local sonic 
values in space and time. This formula can be directly ex- 
tended to the three-dimensional case because of the rotated 
differencing along the streamwise direction: 

A *F = A >< 7 )e. ( 8 ) 

with 



J = A 3 J 



/± r = (^) I>r {l±sign(^, r )cos 3 ^ r } 

± si tr A*(pg)/, r cos 3 0i, r ± pi t r\U m |sin 3 8 - Ajsin 3 0 
pf r = P( t r{l ± sign(l r i, r ) cos 3 0f, r } - A a sin 3 0 

and where Aj = (c m — a^d Aj = ( c m — 

The evaluation of f* r requires many opera- 
tions, but this is not excessively expensive in computations 
because it is a scalar quantity. 


CD Method 

To compare with the upwind formulations, the second- 
order-accurate CD method with the artificial dissipation 
terms 8 is shown here as 
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where k = 0.05 in the present study, 6 2 is a central second- 
difference, cr is a spectral radius of the Jacobian matrix, A, 
and ( ) J+ j = {( )>+i + ( );}/ 2. The parameter a controls 

the strength of the second-order dissipation: 


/3 + 2Af 00 ' , \ 
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\ 4 / 
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(18) 


Note that this model uses the difference of the conservative 
variables directly, while the present upwind formula uses the 
combination of differences, such as total mass flux, pressure, 
and contravariant velocity. 


Results 


Assuming conical flow, the computations become two- 
dimensional. Computations for conical flows are carried out 
in the following manner. The LU-ADI method, 13 which can 
be modified for the conical flow fields, * is used for three 
upwind algorithms and the CD method, which are computed 
explicitly. Laminar flow is also assumed. The third-order 
MUSCL scheme with Koren’s differentiable limiter* is used 
for the two-dimensional computations. 

10° Cone at Zero Incidence 

The first test case is the computation of flow past a 10° 
cone at zero incidence. 3 For a circular cone at zero incidence, 
the solutions become one-dimensional, varying only with the 
angle between conical ray and cone surface. 

Figure 1 shows a comparison of numerical solutions on 
37 grid points. The numerical fluxes are evaluated by three 
first-order-accurate upwind algorithms, ».c., the present for- 
mula (indicated as Goorjian-Obayashi), Roe s FDS and 
Steger- Warming’s FVS formulas, and the CD method. The 
results indicate that the present formula resolves shocks and 
shear layers accurately in the manner of Roe’s formula which 
has a perfect resolution in the one-dimensional case. The 


CD method works fairly well for capturing the boundary 
layer. Steger- Warming’s solution appears most dissipative, 
as pointed out in Ref. 2. 

75° Delta Wing 

The second test case considers a vortical flow field in or- 
der to examine the present formula’s capability for comput- 
ing shear flows. Computations are done for flows past a 
75° delta wing at = 2.8 and a = 16*, for which ex- 
perimental data arc available. 18 Figure 2 shows the model 
geometry, and the typical experimental flow field is shown 
in Fig. 3. For the computations, the conical assumption 
is used. 14 Three grids are used for a grid-refinement study. 
The coarse, medium, and fine grids use 27 (circumferential) 
x 51 (normal to body), 51 x 51 and 99 x 51 points, respec- 
tively, as shown in Fig. 4, where (y,z) corresponds to the 
physical coordinates, not to the conical coordinates (K, Z). 
The vertical solid lines near y = 0.1 and 0.2 in Fig. 4 in- 
dicate the locations of the following comparisons in Figs. 7 
and 8. 

Computations were done with the present method, Roe’s, 
and the CD methods. Figure 5 shows a comparison of den- 
sity contour plots of three numerical solutions on the fine 
grids. Two shock waves can be observed. One is the bow 
shock wave in the windward side of the delta wing, and the 
other is the crossflow shock wave in the leeward side. The 
present and Roe’s methods give similar contour plots for 
those shock waves, but the CD method gives smeared plots. 
(For the CD method, the smoothing coefficient, *, was set 
to 0.1 in the fine-grid case instead of k = 0.05 in the other 
cases because at convergence the solution had numerical os- 
cillations with k — 0.05 in the fine grid.) The low density 
regions at both the primary and secondary vortices also in- 
dicate that both upwind methods give similar solutions, but 
the CD method gives a smeared one. 

A comparison of total pressure contour plots on the fine 
grids is shown in Fig. 6. The primary vortex appears simi- 
larly in both upwind solutions, in respect to its location and 
the contour level, but not in the CD solution. The primary 
vortex appears off from the boundary layer in both upwind 
solutions. But the primary vortex and the boundary layer 
touch each other in the CD solution. The shear-flow region 
separated from the leading edge also shows differences in 
the three solutions. The present formula gives a sharper 
solution than Roe’s in this shear-flow region. The CD solu- 
tion gives the most dissipative solution, again. To confirm 
these observations, solution profiles are compared at y = 0.1 
(approximately on the primary vortex) and y = 0.2 (approx- 
imately on the secondary vortex) as indicated in Fig. 4. 

Figure 7 shows comparisons of total pressure profiles at 
y = 0.1 on the coarse, medium, and fine grids. The large 
gradient near z = 0 corresponds to the boundary layer. In 
the coarse grid case, the CD solution does not have a high 
peak as the other upwind solutions does. This peak becomes 
higher and thus closer to the upwind solutions as the grids 
are refined. The local minimum near z = 0.05 corresponds 
to the primary vortex. The CD method also gives a smeared 
profile. In Fig. 7c, the profile of the CD solution varies con- 
tinuously to the primary vortex. Compared with the CD 
solution, both upwind solutions have a sharper gradient re- 
gion, which confirms the observation in Fig. 6. Roe’s solu- 
tions agree with the present solutions on the medium and 
fine grids but not on the coarse grid. Although not shown in 
this paper, plots of helicity density show that Roe’s solution 
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EXTENSION OF A STREAMWISE UPWIND ALGORITHM 
TO A MOVING GRID SYSTEM 


Shigeru Obayashi,* Peter M. Goorjian, and Guru P. Guruswamy 
Ames Research Center 


SUMMARY 


A new streamwise upwind algorithm has been derived to compute unsteady flow fields with the use 
of a moving-grid system. The temporally nonconservative LU-ADI (lower-upper-factored, alternating- 
direction-implicit) method has been applied for time-marching computations. A comparison of the tempo- 
rally nonconservative method with a time-conservative implicit upwind method indicates that the solutions 
are insensitive to the conservative properties of the implicit solvers when practical time-steps are used. Us- 
ing this new method, computations have been made for an oscillating wing at a transonic Mach number. 
The computed results confirm that the present upwind scheme captures the shock motion better than the 
central-difference scheme based on the Beam-Warming algorithm. The new upwind option of the code 
allows larger time-steps and thus is more efficient, even though it requires slightly more computational 
time per time-step than the central-difference option. 


INTRODUCTION 


A code, ENSAERO, is being developed at Ames using the Euler/Navier- Stokes equations for com- 
puting the unsteady aerodynamics and aeroelasticity of aircraft. The capability of the code has been demon- 
strated by computing vortical and transonic flows over flexible swept wings (refs. 1 and 2). The flow fields 
were calculated by a time-accurate, finite-difference scheme based on central differencing. 

The purpose of this study is to enhance the algorithm capability of the present code. In this re- 
spect, the use of a new upwind scheme in comparison to the current central-difference (CD) scheme is 
investigated. The CD scheme requires an artificial dissipation to stabilize computations. Such artificial- 
dissipation models lead to more dissipative solutions than upwind schemes. In addition, the CD scheme 
is sensitive to the amount of dissipation and needs a specific dissipation coefficient for each case. On the 
other hand, upwind schemes do not require that any coefficient be specified. 

Recently, a streamwise upwind algorithm has been developed and applied to treat steady-state prob- 
lems of transonic flows over wings (ref. 3) and vortical flows over a delta wing (ref. 4) on fixed grids. The 
main feature of the streamwise method is the use of the local stream direction, flow velocity, and pressure 
gradient. The switching of flux evaluations always takes place at sonic values, where shock waves may ex- 
ist. Therefore, this method follows the flow physics more closely than conventional upwind methods based 

* MCAT Institute, San Jose, California. 
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The Cartesian velocity components u, v, and w are nondimensionalized by the free-stream speed of sound 
Coo’, the density p is nondimensionalized by the free-stream density pc*,; and the total energy per unit volume 
e is nondimensionalized by PooC^. The viscous flux vector G v is given by 


with 


G v 


0 

pmiu^ + fm2Cz 
pmi V £ + 3 TU2Cy 
pmi + ^rmCz 

_pmirri 3 + £7712 (Cx u + Cj/ V + C* w ) - 


mi=d +< v 2 + C 2 

m2 = Ci«< + Cv v c + 

m 3 = l -{u 2 + v 2 + w 2 \ + p r(7 "_ j) (c 2 )< 


where Re is the Reynolds number, Pr is the Prandtl number, c is the speed of sound, and J is the trans- 
formation Jacobian. Pressure is related to the conservative flow variables Q, through the equation of state 
for a perfect gas: 


p = (7 - 1) 


e — j(u 2 + v 2 + tu 2 


) 


( 2 ) 


where p is the fluid density and e is total energy per unit of volume of the fluid. 

For the inviscid case, the viscous flux G v is replaced by 0. For the viscous case, the viscosity 
coefficient p in G v is computed as the sum of pi + p t where the laminar viscosity p/ is taken from the 
free-stream laminar viscosity, assumed to be constant for transonic flows, and the turbulent viscosity p* is 
evaluated by the Baldwin-Lomax model. 5 


NUMERICAL ALGORITHM 


The space-discretized form of equation (1) can be written as 


8rQ = 


E <+k 


Ei. 


F,.i 


Ej- 


G 




Ql 


G'i. 


Gl 




A r\ 


AC 


Re 


AC 


(3) 


where a second-order central-difference evaluation is applied to the viscous term. 

The evaluation of the inviscid fluxes is based on the finite-volume cell-centered scheme. To be 
consistent with the finite-difference scheme in ENSAERO, the metrics are defined at each grid point where 
the flow variables are stored. The surface vector of cell-interface, which is necessary for the finite-volume 
formulation, can be obtained by averaging the metrics at the adjoining points. The free-stream preservation 
of this metric evaluation was shown in reference 6. 
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Equation (4) is written in a vector form. If the formula is rewritten in component form, the present 
algorithm can be summarized with a surface vector S and a motion of its centroid, \ t = (x t ,yt, z t ) , as 


F(Qi,Q r ,S 


7+r 


v \ _ gnl / t ?+ , n 1 — \ 


(5) 


where 




ftr 

ft r Hr + k * Pfr 
ft Hr + k V Pfr 

ft + k z P% 

L ftr Hi, r - kt pf r + ( \V m \Ap - V m A 2 )sm 2 6 . 


ft = r [1 ± sign (Vl,r) cos 2 #,, r ] 

± s t ,r A*(pq\ r cos 2 6i, r it pi, r\V m \sin 2 6 - ^-Aisin 2 6 
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where q and V are defined above, and where Ai = (c m — \V m \)^f-, A 2 = (c m — \V m \) p m AV , and 
k t = —k x Xt — k yl/t — k z*t- The averaged state (m) is defined for p, u, v, w, and H by the arithmetic 
average of the left (0 and right (r) states. 

The switches si and s r are defined in the manner of Godunov’s method as follows. For V > 0 , 

5t=i-QC m (6) 

S r = (1 -€m)U -Cr) 

where ^ 

U,m,r = 2 t 1 + Sign ( A/fm.r D] 

and M = q/c. 

A simple way to evaluate the rotation angle is to use cos# = V jq. In supersonic flow fields, 
however, it is important to detect whether the velocity projected to the grid line is beyond the Mach cone. 
Thus, V /q is replaced by M ■ V /q = V /c. If V /c becomes larger than one, cos# is set to one. To avoid 
expansion shocks, the rotation angle is determined by a mixture of averaged (m) and pointwise ( l , r) values: 


cos 2 0i, r = min[ ( 1 —(f)) 


V. 


Vj 2 

+ 

C lr 


1 1 


(7) 
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Approximate Block ADI Method 


The LU-ADI method described in the previous section is nonconservative in time owing to the 
diagonalization. To investigate the significance of this temporal nonconservativeness, an alternative ap- 
proach is considered here. A time-conservative method can be constructed using a block-tridiagonal solver 
similar to the Beam- Warming method. Since true Jacobians of the numerical fluxes of the present upwind 
algorithm are expensive to compute, approximate Jacobians are used here. 


To construct an implicit method for the present upwind algorithm, it is easier to start from the vector 
form (eq. (4)) rather than from the component form (eq. (5)). From equation (4), the first-order-accurate 
flux can be rewritten as 


Fp =- !— ^ { Fj [1 ± sign (V,) cos 2 Oj] ± s ; A*(pg),cos 2 0,e 3; 
2 J j±T 



± P;±^(C;±[ - |V,^l)sin 2 0 e d; . ± . Vj] 


(13) 


To simplify the formula further, the e S;± i. term in the second line of the above formula is replaced by e 3j . 
Then, the Jacobian dF^jdQ can be approximated as 
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' 2 J i± i 1 
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where B = dF/dQ, I is the identity matrix, 
1 dpe. 
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noticed in the unsteady results. The convergence of the unsteady computations to a periodic flow is verified 
by comparing the results between cycles. For all cases considered here, the results for the third cycle give 
pressure profiles that are identical to those of the second-cycle results. Thus, the numerical transient is 
confirmed to disappear within two cycles. 

The present upwind results are also compared with the results obtained by using the existing CD 
method. Reference 13 discusses the CD results obtained from the same test case. The CD method uses 
the diagonal inversions of the Beam-Warming method (ref. 9), and thus it is first-order accurate but non- 
conservative in time. Its artificial dissipation model consists of second- and fourth-order dissipation terms 
controlled by the amount of the second-difference of pressure. 


Inviscid Solutions 


First, the inviscid coarse-grid (91 x 25 x 25 points) solutions are presented. Figure 1 compares the 
computed and measured unsteady, upper-surface pressure coefficients for the real and imaginary parts of 
the first Fourier component. Pressure coefficients are shown for various span locations. The computations 
were done by using the upwind algorithm with the LU-ADI method (UP-LU). The solid, dashed, and dotted 
lines in the figure indicate the results obtained using 1440, 720, and 360 time-steps per cycle of oscillation 
(steps/cycle), respectively. (The time-steps per cycle rate of 1440 corresponds to A t ~ 0 .02 .) The solution 
profiles obtained using 1800 steps/cycle coincided with those using 1440 steps/cycle. Thus, the unsteady 
pressure profiles converged with respect to time-step sizes at 1440 steps/cycle. In addition, because the 
differences between the results with 720 and 1440 steps/cycle are not critical when making comparisons 
with experimental data, the computation with 720 steps/cycle is acceptable for numerical efficiency. 

Figure 2 shows the analogous CD results. In this case, the second- and fourth-order dissipation 
coefficients are fixed at 0.25 and 0.01, respectively (denoted as CD(0.01)). These coefficients are the 
values that were used in references 1, 2, 9 and 13. The CD results show less dependence on time-step 
size than do the UP-LU results, although the profiles are smeared out at the large-gradient region near the 
shock wave and leading-edge regions. In order to check the dependence of the CD method on the amount 
of numerical dissipation, computations were tried with half the amount of dissipation coefficients (0.125, 
0.005). At 360 steps/cycle, the computation became unstable with these coefficients. At 720 and 1440 
steps/cycle, the computations were stable, but the resulting unsteady pressure profiles were still smeared 
out. 


Figure 3 shows the results using the approximate block ADI method applied to the present upwind 
algorithm (UP-BL). The UP-BL method was not stable for computations using 360 and 720 steps/cycle, and 
thus the computations were done using 1080, 2160, and 3600 steps/cycle. The solution profiles obtained 
with 3600 steps/cycle finally showed good agreement with the UP-LU result converged at 1440 steps/cycle. 

Figure 4 compares the three methods, UP-LU, UP-BL, and CD(0.01), fixing the number of time- 
steps per cycle at 1440. Both UP-LU and UP-BL results give similar profiles at the shock motion, although 
the UP-LU method is nonconservative in time. Both upwind methods give better agreement with the 
experimental data in the region of the shock motion than the CD method. 
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shown by oscillations in the pressure profiles of CD(0.01). However, the UP-LU result shows that the UP- 
LU method relaxes this stability requirement. Thus, time-step sizes for the UP-LU method are not limited 
by stability considerations even for the viscous case. 

Figure 7 shows the corresponding results with 1440 steps/cycle. This is the time-step rate sug- 
gested in reference 13. Here, the CD(0.01) solution becomes stable and agrees reasonably well with the 
experimental data. To check whether the amount of dissipation is reasonable, the dissipation coefficients 
were reduced by half, as was done for the inviscid case, but the CD(0.005) computation diverged for the 
viscous case. In figure 7, and analogous to the inviscid results shown in figure 5, the CD solutions approach 
the UP-LU solution as the dissipation coefficients are reduced. Again the UP-LU method is confirmed to 
give a less-dissipative and thus more accurate solution than the CD method. The viscous solutions in fig- 
ure 7 are similar to the inviscid solutions in shown in figure 5 because the experimental flow field does not 
contain strong viscous effect (ref. 12). 

These results also indicate that the CD method is sensitive to both time-step size and to the amount 
of artificial dissipation added. Therefore, users of the CD method are required to find an adequate combina- 
tion of time-step size and dissipation coefficients on a case-by-case basis. In contrast, the UP-LU method 
is robust, and users do not have to find any dissipation coefficient. 

For the viscous case, the UP-LU computation requires 18.9 /isec per grid point per time step at 
a speed of 141 MFLOPS on a CRAY-YMP computer using a single processor, and the CD computation 
requires 17.1 /rsec at a speed of 132 MFLOPS. There is 11% increase of CPU time when the upwind option 
of the code is used. 


CONCLUSIONS 


A new streamwise upwind algorithm has been derived to compute unsteady flow fields with a 
moving grid system. The temporally nonconservative lower-upper-factored, altemating-direction-implicit 
(LU-ADI) method has been applied to compute flow over an oscillating wing at a transonic Mach number. 
A comparison of the temporally nonconservative method with a time-conservative version of the upwind 
scheme indicates that the solutions are insensitive to the time-conservativeness of the implicit solvers when 
practical time-step sizes are used. The temporally nonconservative upwind method was found to be more 
stable and twice as efficient as the time-conservative block ADI scheme, even though the block ADI scheme 
employed an approximate form of the Jacobians. Comparisons have been made between results obtained 
with both upwind schemes, with experimental measurements, and with computed results obtained using 
the existing central-difference method. 

Comparisons with experimental data show that the upwind algorithm predicts the shock motion 
better than the central-difference (CD) method. The CD solutions are also found to be sensitive to the 
amount of numerical dissipation, and as a result, the dissipation coefficients must be specified case by case. 
In comparison, the present upwind method does not require a dissipation coefficient. Thus, the method that 
combines the streamwise upwind and LU-ADI methods is proposed for practical computations. This new 
upwind method is incorporated in the aeroelastic code ENSAERO. The new upwind option of the code 
allows larger time-steps and thus is more efficient, even though it requires slightly more computational 
time per time-step than the CD option. 
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Abstract 

A new streamwise upwind algorithm has been derived 
to compute unsteady flows with a moving grid system and 
applied to compute flows over oscillating wings at tran- 
sonic Mach numbers. Comparisons have been made be- 
tween results obtained from this upwind algorithm, using 
both temporally nonconservative- and conservative-implicit 
methods, with the results obtained from a central-difference 
method, and also with experimental data. The results show 
(1) the efficiency and practicality of the temporally noncon- 
servative implicit solver and (2) the robustness and accuracy 
of the upwind method for unsteady computations compared 
to the central-difference method. 


Introduction 

In the last two decades, there have been extensive de- 
velopments in computational aerodynamics, which consti- 
tutes a major part of the general area of computational fluid 
dynamics. Such developments are essential to advance the 
understanding of the physics of complex flows, to comple- 
ment expensive wind-tunnel tests, and to reduce the overall 
design cost of an aircraft, particularly in the area of aeroe- 
lasticity. 

Aeroelasticity plays an important role in the design and 
development of aircraft, particularly modem aircraft, which 
tend to be more flexible. Several phenomena that can be 
dangerous and limit the performance of an aircraft occur be- 
cause of the interaction of the flow with flexible components. 
For example, an aircraft with highly swept wings may ex- 
perience vortex-induced aeroelastic oscillations. 1 Also, un- 
desirable aeroelastic phenomena due to the presence and 
movement of shock waves occur in the transonic range. 
Aeroelastically critical phenomena, such as a low transonic 
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flutter speed, have been known to occur through limited 
wind-tunnel tests and flight tests. 

Aeroelastic tests require extensive cost and risk. An 
aeroelastic wind-tunnel experiment is an order of magnitude 
more expensive than a parallel experiment involving only 
aerodynamics. By complementing the wind-tunnel exper- 
iments with numerical simulations, the overall cost of the 
development of aircraft can be considerably reduced. In or- 
der to accurately compute aeroelastic phenomenon it is nec- 
essary to solve the unsteady Euler/Navier-Stokes equations 
simultaneously with the structural equations of motion. 

At Ames a code, ENSAERO, is being developed for 
computing the unsteady aerodynamics and aeroelasticity of 
aircraft and it solves the Euler/Navier-Stokes equations. The 
capability of the code has been demonstrated by computing 
vortical and transonic flows over flexible swept wings. 2,3 
The flow fields were calculated by a time-accurate, finite- 
difference scheme based on central differencing. 

The motivation of this study is to enhance the algorithm 
capability of the present code. Toward this goal, the use of 
a new upwind scheme in comparison to the current central- 
difference (CD) scheme is investigated in this paper. The 
CD scheme requires artificial dissipation to stabilize com- 
putations. In general, such artificial dissipation models lead 
to more dissipative, and thus less accurate, solutions than 
upwind schemes. In addition, the CD scheme is sensitive 
to the amount of dissipation and it is necessary to specify a 
dissipation coefficient on a case by case basis. On the other 
hand, upwind schemes do not require any coefficient to be 
specified. 

Among upwind algorithms, a streamwise upwind al- 
gorithm has recently been developed and applied to steady- 
state problems of transonic flows over wings 4 and vortical 
flows over a delta wing 5 on fixed grids. Most multidimen- 
sional upwind algorithms are first constructed in one dimen- 
sion and then extended to multidimensions by applying the 
one-dimensional procedure in each coordinate direction. By 
comparison, the present method uses the local stream direc- 
tion, flow velocity, and pressure gradient to construct the 
upwinding. The switching of flux evaluations always takes 
place at sonic values, where transonic shock waves may be 
located. Therefore, this method follows the flow physics 
more closely than the coordinate upwind methods. The com- 
puted results confirmed the higher resolution of the present 
algorithm over the CD method as well as over other upwind 
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For a moving system. 


Vmot/in* = Vt + T) X U + TfoV + 7],^ 

= T? x (u - X t ) + 77j,(v - y t ) + 17 jf( xtf - 2t) 
= Vtj $ 


where = — rj x x t — r^yt — rjt z t * Then the formulas for a 
fixed system can be rewritten for a moving system. Note that 
the present algorithm uses x t » Vt> and zt to obtain the flow 
velocity for upwinding in the streamwise direction, instead 
of fjt. and Ct to compute the contravariant velocity for 
upwinding in the coordinate direction. 

The present algorithm can be summarized with a sur- 
face vector S = ^ and a motion of its centroid, = 


F(Q«,Q r ,S, + *,x v+ p = !^!( jl K + Ft 1 " M3oo) 

(4) 

where 
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A, =(a n -\V n \ 

°m 

A 2 =(o m -|V' m |)AnAV 


with A • = - T - -i, V = + k,w, kt = - 

k x = r),/|Vtj|, and so on. The subscripts l, m 
and r denote the left, averaged and right state of the flow 
variables. The averaged state is defined for p, u, v, w, and 
H by the arithmetic average of the left and right states. The 
basic scheme is first-order accurate with l = j and r = j + 1. 
The term, JbtQoo. subtracts the free-stream for time-metrics. 


The switches si and s T are defined in the manner of Go- 
dunov’s method as follows. For V > 0 , 


where 


Sl — 1 £l Em 

Sr = ( I C m ) ( 1 £r) 

£t.m.r = j{ 1 + sign( Mf mj - 1)} 


(5) 


and A f = q/a. 

A simple way to evaluate the rotation angle is to use 
cos0 = V /q. In supersonic flow fields, however, it is im- 
portant to detect whether the velocity projected to the grid 
line is beyond the Mach cone. Thus, V/q is replaced by 
M -V/q = V/a. If V /a becomes larger than one, cos0 is 
set to one. To avoid expansion shocks, 5 the rotation angle 
is determined by a mixture of averaged (m) and pointwise 
(l t r) values: 

„ y 2 V 2 

cos 2 0ir = min[ (1 - 4>)~r + 1 1 ( 6 ) 

a m *l t r 

The following relation is used for evaluating <f> e (0, 1] in 
this paper because of the smoothness: 

*=max[ J^L(i_J_ {7 _i + (7+ i ) W}) i 0] (7) 

where pi and pi denote upstream and downstream pressures, 
respectively. The sine is determined by an arithmetic aver- 
age of the cosines: sin 2 0 = 1 - y(cos 2 0i + cos 2 0 r ). 

Higher-order schemes are constructed from a one- 
parameter family, k , of interpolations of the primitive vari- 
ables, p, u, v, w y and p. For example, 


Pi = {l+^[(l -*)V+(l + iOA]} P> 

p r = { i _^l[( 1 + | t )V+( 1 -«)A]>p /+1 (8) 


where V and A are backward and forward difference op- 
erators, respectively. 9 For the third-order scheme, k - $•, 
Koren’s differentiable limiter 10 is used in this paper. The 
limiter y> is calculated as 


. _ 3Vp, A p ; + e 

' 2(A Pj - Vp ; ) 2 + 3Vp ; Ap ; + £ 


(9) 


where a small constant £, e = 10 ~ 6 typically, is added to 
prevent the division by zero. The same formulas are used 
for the other primitive variables. 


LU-ADI Method 

The time marching method used for the present upwind 
scheme is the LU-ADI factorization method proposed by one 
of the present authors. 7 The LU-ADI method is a compro- 
mise of ADI and LU factorization. This method applied to 
Eq. (3) is written as 

(T ( L A D A [/ A Tf') (T'LbDbUbT- 1 ) 
ftLcDcUcTf 1 ) = AtR' ( 10) 


794 



Unsteady computations are started from the corre- 
sponding steady-state solution. The convergence of the un- 
steady computations to a periodic flow is verified by com- 
paring the results between cycles. For all cases presented 
here, the third-cycle results give identical pressure profiles to 
those of the second-cycle results. Thus the numerical tran- 
sient is confirmed to disappear within two cycles. 

First, the LU- ADI upwind method was applied to com- 
pute this inviscid flow on the coarse (91 x 25 x 25 point) 
grid using 360, 720, 1080, 1440, and 1800 time steps per 
cycle of oscillation (steps/cycle). The unsteady pressure 
profiles on the wing surface that were obtained using 1800 
steps/cycle coincided with those using 1440 steps/cycle. The 
computation converged with respect to time-step sizes at 
1440 steps/cycle. Except for the result using 360 steps/cycle, 
differences between the other results were small. Thus, 720 
steps/cycle is acceptable for practical computations. 

The CD method showed less dependence on time-step 
sizes than the LU-ADI upwind method. The CD computa- 
tion converged at 720 steps/cycle. The block-ADI upwind 
method showed more dependence on time-step size. The 
solution profiles obtained with 3600 steps/cycle finally gave 
a good agreement with the LU-ADI upwind result obtained 
at 1440 steps/cycle. The block-ADI computations were not 
stable for 360 and 720 steps/cycle. 

Figure 1 shows the comparison of real and imaginary 
parts of the first Fourier component between the computed 
and measured unsteady upper-surface pressure coefficients 
of the wing at various spanwise locations, using the three 
methods with 1440 steps/cycle. Both upwind results give 
similar profiles of the shock motion, although the LU-ADI 
upwind method is nonconservative in time. Both upwind 
results give crisper profiles (narrower peaks) at the region of 
the shock motion than the CD result. 

These results indicate that the temporally nonconserva- 
tive LU-ADI method can be used for unsteady computations 
even with moving shock waves when a large enough number 
of time steps (i.e., small A t) is used per cycle. In the present 
case, any number greater than 720 will give practical results. 
The CD method also uses the temporally nonconservative 
diagonal form. However, the result differs from both up- 
wind results. This indicates that the solution depends on the 
numerical dissipation more than the time-conservative prop- 
erties of the methods. The block-ADI method requires twice 
as much CPU time as the LU-ADI method, but its accuracy 
does not appear to compensate for the increased computa- 
tional time. Thus, the block-ADI method will be dropped in 
the following computations. 

Next, the inviscid computations were repeated on a 
finer (151 x 25 x 34 point) grid to check the grid depen- 
dency. Figure 2 shows the comparison of unsteady pressures 
using the LU-ADI upwind and CD methods. Both compu- 
tations use 720 steps/cycle, as suggested in the coarse-grid 
case. Compared with the coarse-grid solution in Fig. 1, the 
CD solution tends to converge to the upwind solution due 
to the grid refinement, although the profiles are still slightly 


smeared at the peak. Therefore, the upwind method is con- 
firmed to give a less dissipative and thus more accurate so- 
lution than the CD method. 

F-5 Wing 

The second test considers unsteady viscous flows over 
an F-5 wing which has an aspect ratio of 2.98, a taper ratio 
of 0.31 and a leading edge sweep angle of 31.92°. Com- 
putations were made using two grids: the coarse and fine 
grids containing 121 x 25 x 25 points and 151 x 25 x 
30 points, respectively. Figure 3 shows the F-5 wing and 
the grid distributions at the root section of the fine grid. Fig- 
ure 4 illustrates the motion used in the experiment conducted 
at the National Aerospace Laboratory of the Netherlands. 16 
The wing is pitching about an axis located at the 50% root 
chord, and the pitching axis is normal to the wing root Ref- 
erence 3 discusses the CD results applied to the same test 
case. 

The test cases are chosen at Mo© = 0.896 and 1.328 
where the measured steady and unsteady data are given in 
Ref. 16. All F-5 wing cases are computed at Re = 9 x 10 6 
based on the root chord. At these flow conditions, the coarse- 
grid solution gave y* < 13.3 at the first shell of points above 
the wing surface. For both the steady and unsteady cases, the 
mean angle of attack a m is 0°. The unsteady flows are com- 
puted at a reduced frequency k = 0 .550 and a pitch ampli- 
tude a = 0.11° for Moo = 0.8%. Similarly,k= 0.3% and 
a = 0.22° for Moo = 1 328. The Baldwin-Lomax eddy- 
viscosity model is used to compute the turbulent viscosity 
coefficient 

The time-step dependency of the upwind method was 
checked on both the coarse and fine grids at Mo© = 0 .8% . 
The unsteady pressure profiles converged with respect to 
time-step sizes at 1800 steps/cycle for the coarse grid and 
2160 steps/cycle for the fine grid. For computational effi- 
ciency, 1440 and 1800 steps/cycle are deemed acceptable 
for the coarse and fine grids, respectively. All results shown 
here were computed with 1800 steps/cycle. 

Moo = 0.8% 

Figures 5 and 6 show the comparisons of the computed 
steady and unsteady pressures with the experimental data at 
Mo© * 0 .8% using the upwind and CD methods on the 
coarse (121 x 25 x 25 point) grid. For the steady case, 
both numerical solutions show shock-free profiles in agree- 
ment with the measured data (note that CJJ = —0 1% for 
Moo = 0 .8% ). For the oscillatory case, unsteady pressure 
peaks of the upwind solution indicates that the unsteady mo- 
tion of the wing produces a shock wave. The measured data 
also show the peaks except at the 20% semispan location. In 
contrast, the unsteady pressure peaks of the CD solution are 
smeared out and the resulting profiles agree less favorably 
with the experiment except at the 20% semispan location. 
The computation assumes symmetry at the wing root instead 
of a side wall which may affect the experimental data at the 
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AGARD rectangular wing 
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□ Exp, ref. 15 
94% semispan 



77% semispan 



50% semispan 




F| g. 2 Comparison of computed inviscid upper surface unsteady pressures between the LU-ADI upwind and central- 
dilference method over the rectangular wing with the fine grid. 
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Steady 








F-5 wing 

= 0.896 

a - 0° (pitching down) 
k = 0.550 
Re = 9 X 10 6 
Grid: 151 X 25 X 30 
70% semispan 



b) 


Fig. 9 Comparison of computed instantaneous pressure contours between the upwind and central-difference methods at 
Mao = 0 .896 over the F-5 wing with the fine grid, a) Upwind result, b) central-difference result 
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Introduction 

Body-conforming coordinates transformation of a fluid conservation- law form is 
generally used in computational fluid dynamics. The metrics associated with the coordi- 
nates transformation are required to satisfy certain geometric identities to maintain the free 
stream. 1 These metrics are called free-stream capturing (or preserving) metrics. So far, nu- 
merical techniques are known to capture the free-stream on stationary grids. 2-4 However, 
the extension of the free-stream capturing metrics to moving grids is not straightforward. 
The error introduced by the time metrics has been overlooked because it is negligible in 
most cases, but it can be significant in certain applications such as helicopter rotor flow 
fields. 5 

Rigorous formulations to avoid this error were suggested in Ref. 1, and demon- 
strated, for example, in Ref. 6. Based on the work in Ref. 1, the present study describes 
detailed formulas for constructing the free-stream capturing metrics in space and time on 
both the finite- volume (FV) and finite-difference (FD) framework. The error introduced 
by the inconsistent time-metric term is also evaluated. 


* Research Scientist, Applied Computational Fluids Branch, MCAT Institute, San 
Jose, California 95127. 
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Using Eqs. (4) and (6), one obtains 



r x ndS = 0 


(7) 


The conservation of volume for a time-varying cell is given by 


V(t 2 )-V(t 0 = 



n • v c dS dt 


(8) 


The free-stream preservation due to these geometric identities can be demonstrated by 
substituting Qoo, Fstoc and Eqs. (2) to (4) to Eq. (1): 


Qoo[V(h) - V(ti)] = - /[ F sioo • ( / ndS) - <?oo(v o - ft X ro) ■ ( / n dS) 

Jt ! Js Js (g) 

- Qo oft • (j> r x ndS) - Qoo(j> n • v c dS)}dt 


The geometric identity, Eq. (5), suffices to capture the free stream in a fixed coordinate 
system, where most of steady-state computations are carried out, and in a moving coor- 
dinate system without rotation (ft = 0). When the grid is moving with rotation (ft ^ 0), 
the second geometric identity, Eq. (7), is to be satisfied. For the general motion of grid 
with changing cell volume, the third geometric identity, Eq. (8), is also required. 

The geometric identity, Eq. (5), preserves the free stream at any instance t when 
v = 0. Thus, the time-differential form of Eq. (1) is often used in the FV formulation. 
However, if the grid moves, the geometric identities have to be satisfied correctly in the 
integral form. 


Free-Stream Capturing in the Inertial Frame 

To preserve the free-stream perfectly with a moving grid, Ref. 1 suggests to consider 
the rigorous FV formulations in space and time. One of the rigorous FV formulations with 
a grid velocity expressed in the inertial frame is shown here in detail. 

Assuming that v is given and then by substituting Q oo and F,^ to Eq. (1), one 
obtains 


Qco[V(h) — V(fl)] = —Fstoo 



n dSdt + Q oo 



n • vdSdt 


( 10 ) 
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where Sn* 5<5 and Si22'i' are the surface vectors in space and time domain. Note that this 
formula requires only the difference of the positions of grid points between 1 i and £ 2 , not 
the grid velocity v itself. Then, instead of Eqs. (6) and (8), one obtains 

A V = Y, y s (18) 

cell 

This identity does not mean to satisfy Eq. (7) but does satisfy Eq. (6) in the time- 
integral form. Therefore, Eq. (18) leads to the perfect free-stream capturing with the 
use of Eq. (11). 

Free-Stream Capturing in the Non-Inertial Frame 

It is convenient to use the non-inertial frame for certain applications. Thus, the 
FV formulation with a grid velocity expressed in the non-inertial frame is shown next. 
The analysis also provides a deep insight for the free-stream capturing because it considers 
three types of motion given in Eqs. (3) and (4) separately. The discretized forms of the 
geometric identities in the FV method can be expressed as 


E 5n = 0 

(19) 

cell 

Sr x n = 0 

(20) 

cell 

av = J2 v sc 

(21) 


celt 


where Vs c is obtained from Eq. (16) by replacing v with v c (see Ref. 1 for more details). 

Reference 1 introduces the axea moment M = J s r x ndS so as to satisfy the 
discretized geometric identity, Eq. (7), on the hexahedron: 

Mi 562 = 1*165 X Si 65 + 1*126 X S126 (22) 

where ri 65 = •5(1*1 + r6 + rs), S165 = 5 ( 1*6 — ri) x (rs — ri), and so on. Note that 
M 1562 ^ 1*1562 x Si 562 - The expression, ri 562 x S 1562 , is not well-defined for computing 
area moment. In contrast, Eq. (22) is well-defined. To see these, let r<j>, r^, and Sa 
be ri562, 1*165, |Si 562 |n and |Si6s|n, respectively. After simple algebraic manipulation, one 


5 



The other conceptually different way is to use S a always; that is, to regard the 
hexahedral cell as dodecahedron or to divide the hexahedral cell into tetrahedron. Then, 
instead of Eq. (14), one will obtain 

J>a = 0 (28) 

cell 

for either dodecahedron or tetrahedron in addition to Eq. (24). The resulting metric terms 
will preserve the free stream. The use of the tetrahedral cell allows the most compact and 
consistent metric formulation. Note, however, that the use of the tetrahedron results in 
unstructured-grid formulations. 

Finite-Difference Formulation 
Geometric Identities in the Finite-Difference Formulation 

The analysis of the FD method can be simplified with the aid of the above discussion 
of the FV method. The FD formulation has to be derived from the integral form, Eq. (1). 
Again from Ref. 1, the differential form for Eq. (1) can be written with a generalized 
coordinate transformation, 

r = r(^,r;,C,r), t - r (29) 

ELS follows: 

Qr + + Frj + = 0 (30) 

where subscripts indicate partial differentiation, 

Q = QV, E = S* • F, F = S" • F, G = S c • F 

and where 

S ^ = r^ x r c , S'* = r c x r^, = r^ x r v (31) 

and 

V = r € • r, x r^- (32) 
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where 8 denotes the difference operator. After scaling by one-fourth to adjust the area from 
double-sized to regular cell, the first component of the above expression can be written as 

y = (^<A2 /)(/V c z) - (/icMX^Ay) ( 38 ) 

where denotes the arithmetic averaging operator. Thus, Eq. (12) is equivalent to the 
consistently differenced metrics in Ref. 3 that are based on the averaging procedure so as 
to satisfy the differential chain rules numerically. 

The main discrepancy between the FV and FD formulations appears in the defini- 
tion of cell volume. The cell volume defined by Eq. (13) is different from the one defined 
by the discretized form of Eq. (32) because the FD method does not use the cell concept. 
Nevertheless, Eq. (13) can be applied to the FD method with a scaling factor of one-eighth 
instead of Eq. (32). Then, the FV space metrics on the double-sized cell become identical 
to the FD ones. 

The FD time-metric evaluation is also considered from the FV point of view. It is 
easily found that the time-metric evaluation, Eq. (34), will not maintain the free stream 
even with the use of the free-stream capturing metrics, Eqs. (11) or (12), in case of a 
rotating frame because of Eq. (23). Also, it can be shown that such inconsistent time 
metrics do not satisfy GCL. The discretized form of GCL can be written as 

A V = Arfo(S« • r r ) + ^(S* • r r ) + *<( S c • r r )] (39) 

Let the grid move in the rigid rotation, that is, V T = 0 and r r = (lxr. Then the left-hand 
side of Eq. (39) is zero. But the right-hand side results in ^(rxS^)+^,(rxS' ? )-|-^(rxS < *) ^ 
0 (Eq. (23) appears again). This indicates that the use of the GCL condition, Eq. (39), 
for computing AV can be erroneous. In other words, the GCL condition, Eq. (39), is 
necessary to preserve the free stream, but not sufficient to construct consistent metrics in 
space and time. 

It is easily found that the time integration of Eq. (39) from i\ to t 2 results in 
Eq. (18) and thus both equations are equivalent. Therefore, the consistent time metrics 
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Abstract 

Unsteady Navier-Stokes computations have been 
carried out for simulating transonic flows over a clipped 
delta wing undergoing oscillatory and ramp motions, 
including flexibility. The implicit upwind algorithm has 
been validated by comparing the solutions with exper- 
imental data for the oscillatory pitching motion cases. 
The numerical and experimental results agree well at 
moderate angles of attack, where a leading-edge vortex 
develops. The ramp motion cases have demonstrated 
the effects of unsteadiness of the flow field and struc- 
tural flexibility on the wing responses. For the 10° 
ramp motion, a vortex breakdown is observed. The 
interaction with the shock wave plays an essential role 
in the process of the breakdown observed in the present 
calculation. 


Introduction 

In the last decade, there have been extensive de- 
velopments in computational fluid dynamics toward 
the prediction of the three-dimensional flow field about 
complex geometries at moderate and high angles of at- 
tack. Among the characteristics of flows over aircraft, 
the behavior of the flow over delta wings is of strong 
interest for high-speed aircraft because of the nonlin- 
ear lift increase due to the leading-edge vortex. Effects 


^Research Scientist, MCAT Institute, member AIAA. 
**Research Scientist, AIAA Associate Fellow. 
Copyright © 1991 by the American Institute of Aero- 
nautics and Astronautics, Inc, No copyright is asserted 
in the United States under Title 17, U.S. Code. The 
U.S. Government has a royalty-free license to exercise 
all rights under the copyright claimed herein for Gov- 
ernmental purposes. All other rights are reserved by 
the copyright owner. 


of flexibility can further influence the nature of flows 
on such wings. Steady-state flow problems associated 
with delta wings have been widely investigated compu- 
tationally (for example, Refs. 1 to 3). Several advanced 
studies have also been performed for unsteady vortical 
calculations at subsonic and supersonic Mach numbers. 
For example, Ref. 4 presented a conical flow computa- 
tion on a rigid wing and Ref. 5 presented a subsonic 
three-dimensional computation on a flexible wing, both 
involving vortical flows. 

Numerical methods can play an important role in 
complementing expensive wind tunnel tests, particu- 
larly in the area of aeroelasticity. An aeroelastic wind 
tunnel experiment is an order of magnitude more ex- 
pensive than a similar rigid-body experiment involv- 
ing only aerodynamics. By complementing the experi- 
ments with numerical simulations, the overall cost of 
the development of aircraft can be considerably re- 
duced. Thus development of a numerical method is 
desired for simulating aeroelastic phenomenon. To ac- 
curately compute aeroelastic phenomenon it is neces- 
sary to solve the unsteady Euler/Navier-Stokes equa- 
tions simultaneously with the structural equations of 
motion. 

Recently a code, ENSAERO, was developed to 
compute aeroelastic responses by simultaneously in- 
tegrating the Euler/Navier-Stokes equations and the 
modal structural equations of motion using aeroelasti- 
cally adaptive dynamic grids. 6,7 An upwind algorithm 
was implemented into the code and the resulting code 
was successfully applied to compute transonic flow’s 
over a typical fighter-type wing undergoing oscillatory 
motion. 8 

The purpose of this paper is to examine the ca- 
pability of the present numerical method by simulating 
unsteady transonic flows on a clipped delta wing, which 
contain a leading-edge separation. So far, unsteady 
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physics more closely than the coordinate upwind meth- 
ods. The computed results confirmed the higher resolu- 
tion of the present algorithm over the central -difference 
method as well as over other upwind methods . 13 In this 
work, the streamwise upwind algorithm is applied to 
compute the inviscid cell- interface fluxes. A second- 
order central-difference evaluation is applied to the vis- 
cous term. The complete algorithm can be found in 
Refs. 8 and 13. 

An implicit method is used for the time integra- 
tion because the computational efficiency of the method 
is critical for expensive unsteady viscous calculations. 
The method chosen here is the LU-ADI (lower-upper 
factored, alternating direction implicit) method 14 be- 
cause it requires only scalar bidiagonal matrix inver- 
sions. See Refs. 15 and 16 for additional details. 


Aeroelastic Equations of Motion 

The governing aeroelastic equations of motion of a 
flexible wing axe solved using the Rayleigh-Ritz 
method. In this method, the resulting aeroelastic dis- 
placements at any time are expressed as a function of a 
finite set of assumed modes. The contribution of each 
assumed mode to the total motion is derived by La- 
grange J s equation. Furthermore, it is assumed that the 
deformation of the continuous wing structure can be 
represented by deflections at a set of discrete points. 
This assumption facilitates the use of discrete struc- 
tural data, such as the modal vector, the modal stiff- 
ness matrix, and the modal mass matrix. These can 
be generated from a finite-element analysis or from ex- 
perimental influence-coefficient measurements. In this 
study, the finite-element method is used to obtain the 
modal data. 

It is assumed that the deformed shape of the wing 
can be represented by a set of discrete displacements at 
selected nodes. From the modal analysis, the displace- 
ment vector {d} can be expressed as 

w = [4>m ( 4 ) 

where [ 0 ] is the modal matrix and {g} is the general- 
ized displacement vector. The final matrix form of the 
aeroelastic equations of motion is 

[MM + [G]{g} + [*]{?} = {F} (5) 

where [M], [G], and [K] are modal mass, damping, and 
stiffness matrices, respectively. {F} is the aerodynamic 
force vector defined as (^)pf^^[0] r {A]{ACp} and [A] 
is the diagonal area matrix of the aerodynamic control 
points. 


The aeroelastic equation of motion (Eq. (5)) is 
solved by a numerical integration technique based on 
the linear acceleration method . 17 

Aeroelastic Configuration Adaptive Grids 

One of the major deficiencies in computational 
aerodynamics using the Navier- Stokes equations lies 
in the area of grid generation. For steady flows, the ad- 
vanced techniques such as zonal grids 18 are being used. 
Grid generation techniques for aeroelastic calculations, 
which involve moving components, are in early stages of 
development. In Ref. 7, aeroelastic configuration adap- 
tive dynamic grids were successfully used for computing 
time-accurate aeroelastic responses of swept wings. In 
this work, a similar technique is used. 

Results 

Numerical schemes used for flow calculations in 
aeroelasticity must guarantee the correct calculation of 
amplitude and phase of unsteady pressures. To verify 
the accuracy of the present code for simulating the com- 
plicated flow field containing a leading-edge vortex and 
a shock wave, test cases are chosen from the experiment 
on a clipped delta wing undergoing prescribed pitching 
motion . 9 Since the experiment was conducted using a 
Freon test medium, the ratio of specific heats 7 is set 
to 1.135 in the following computations. 

Steady Pressures 

Steady-state calculations have been performed to 
examine the validity of the numerical procedure and 
the computational grid. The model planform geome- 
try is shown in Fig. 1 . The wing has a leading-edge 
sweep angle of 50.4° and a 6 %-thick circular-arc airfoil 
section. The lines, A, B, and C, indicate the spanwise 
locations of the pressure orifices in the experiment. Fig- 
ure 2 shows the grid generated algebraically in the C-H 
topology. The f , r/, and £ coordinates represent the 
chordwise, spanwise, and normal (to the wing surface) 
directions, respectively. The grid contains 151 x 25 x 34 
points. 

The present grid is fairly coarse compared with the 
grids used in the typical steady-state computations . 1 " 3 
The number of grid points was determined to compro- 
mise the accuracy and the total computational time. 
Unsteady computations require more computational 
time than the steady-state computations. With the 
present grid, a typical unsteady case can be computed 
within 5 hr by using a Cray YMP computer with a 
single processor. (The code requires about 19 /isec per 
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at the 68% section in Fig. 8 shows no effect due to the 
leading-edge vortex, it is inconsistent with the other 
data. Overall, the numerical results show fairly good 
agreement with the experimental data. Again, the dif- 
ference of shock motion due to the Mach numbers and 
the difference of the vortex motion due to the angles of 
attack are successfully captured. 

Figure 12 shows the comparison of unsteady pres- 
sures between the modified and unmodified Baldwin- 
Lomax turbulence models for the fourth case (Fig. 11). 
To see the difference clearly, the unsteady pressures are 
plotted in the magnitude and phase angle here. The 
improvements due to the modification of the turbu- 
lence model are seen in both magnitude and phase angle 
where the leading-edge vortex exists. The results also 
confirm that the present grid has reasonable resolution 
for the present unsteady flow fields, even though the 
grid is fairly coarse. Further improvements and grid 
refinements will be needed to achieve a better agree- 
ment at the outboard sections. 

Rigid and Flexible Ramp Motions 

In maneuvering, aircraft often undergo rapid ramp 
motions. During such motions, flow unsteadiness and 
wing flexibility play important roles. In this section, the 
applicability of the present development to computing 
such flow fields is demonstrated. 

Computations axe performed for rigid and flexible 
wings in ramp motion. Structural properties of the 
wing were selected to represent a typical fighter wing. 
Figure 13 shows the mode shapes and the frequencies 
of the first four normal modes for the clipped delta 
wing used in the following computations. The dynamic 
pressure is set to be 1.0 psi . Test cases consider 4 and 
10° ramp motions for both rigid and flexible wings. 

4° Ramp Motion 

Figure 14 shows the comparisons of the sectional 
lift responses between the rigid and flexible wings at 
Moo = 0.9 and Re c = 15 x 10 6 for the wing ramping 
up to 4° with a pitch rate of A = 0.04. The pitch rate 
A is defined as qc/Uoq . The data are plotted at 34%, 
54%, 68%, and 90% semispan sections. The unsteady 
computations are started from the converged steady- 
state solution at 0° angle of attack. 

In the rigid-wing case, the lift responses at the in- 
board sections settle down quickly after the ramp mo- 
tion stops, and the flow approaches the steady-state 
values. At the 90% section, however, the lift contin- 
ues to increase for a short period after the ramp mo- 
tion stops. This corresponds to the movement of the 
leading-edge vortex as indicated in the corresponding 


pressure history shown in Fig. 15. In Fig. 15. the 
pressure distributions are plotted every 100 time steps. 
The ramp motion ends at 600 time steps (t ^ 0.04 sec), 
which corresponds to the sixth plot from the bottom. 
At the 90% section in Fig. 15d, a leading-edge vortex is 
formed and lifts off from the wing surface. This results 
in the increase and then a decrease of the sectional lift 
shown in Fig. 14. At the other inboard sections in Fig. 
15, the pressure distributions become similar to those 
in Fig. 6. i.e., the leading-edge vortex remains without 
interacting with the shock wave. Thus, the correspond- 
ing sectional lift stays nearly steady with minor fluctu- 
ations due to the movement of the vortex, as seen in 
Fig. 15a. 

In contrast to the rigid-wing case, Fig. 14 shows 
that the sectional lift responses of the flexible wing are 
oscillatory. The primary effect of the flexibility is the 
reduction of the lift at most of the spanwise sections 
due to the reduced angles of attack. Figure 16 shows 
the corresponding pressure history plots. The plots in 
Fig. 16c and d show the shedding of the leading-edge 
vortex in comparison to the plots of the rigid wing in 
Fig. 15c and d. In addition, Fig. 14 shows that the 
90% sectional lift response of the flexible wing has the 
low-frequency primary oscillation perturbed by a high- 
frequency vortex shedding. The low-frequency primary 
oscillation corresponds to the structural oscillation as 
showm later. The high-frequency perturbation corre- 
sponds to the shedding of the vortex that can be seen 
in the pressure plots of Fig. 16d. 

Computations (not shown) were also carried out 
for the flexible-wing ramping from 0° to 3° and from 
0° to 5° at the same pitch rate as the 4° case. The 4° 
case show's the largest high-frequency perturbation. In 
the 3° case, the vortex is not strong enough to disturb 
the lift response. The vortex shedding is found only at 
the 90% section. In the 5° case, the vortex lifts off from 
the wing surface so that the structural oscillation does 
not cause the perturbation seen at the lower angles of 
attack. The 10° case discussed in the next section does 
not show the perturbation either. 

An aeroelastic reduction of the local angle of attack 
results in a delay of the lift increase of the flexible wring 
for a short period after the ramp motion stops (about 
0.04 < t < 0.11 in Fig. 14). Figure 17 shows the in- 
stantaneous density contour plots of the 90% spanwise 
section at 800 time steps (t ~ 0.05 sec). The contours 
are plotted at intervals of 0.02 (bold line indicates free- 
stream density, 1.0). The plots show that the leading- 
edge vortex on the flexible wing is w r eaker than that 
on the rigid wing. Figure 18 shows the corresponding 
deformation of the flexible wing. The solid lines indi- 
cate the corresponding surface of the rigid wing. The 
actual displacement of the leading edge at the wing tip 
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Concluding Remarks 

In this paper, a computational procedure for com- 
puting the unsteady transonic flows associated with the 
leading-edge vortex on a clipped delta wing, includ- 
ing flexibility, has been presented. The procedure is 
based on a time-accurate computational method com- 
bined with the use of aeroelastically adaptive dynamic 
grids. The flow is modeled using the Navier-Stokes 
equations. The flow equations are coupled with the 
structural equations to account for the flexibility. The 
numerical procedure has been verified through the com- 
parisons with the experiment for the unsteady pitching 
cases on the rigid clipped delta wing. The main flow 
structures are successfully captured using a fairly coarse 
grid. 

The ramp motion cases have demonstrated the ef- 
fects of unsteadiness of the flow field and flexibility of 
the wing. The primary effect of the flexibility is the re- 
duction of the lift due to the deformation of the wing. 
Interaction of the leading-edge vortex with the shock 
wave has significant effects on the wing responses. For 
the 4° ramp motion, the vortex shedding occurs at the 
wing tip due to the flexibility. For the 10° ramp motion, 
a possible vortex breakdown is observed. The interac- 
tion with the shock wave plays an essential role in the 
process of the breakdown observed in the present calcu- 
lation. A further grid-refinement study will assess the 
validity of the present observations. 
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Fig. 3 Comparison of steady pressures with experiment, = 0.88, a = 3.04, Re c = 17.2 x 10 6 . 
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Fig. 4 Comparison of steady pressures with experiment, = 0.88, a = 4.03, Re^ = 17.3 x 10®. 
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Fig. 7 Comparison of steady pressures between the standard and modified Baldwin-Lomax turbulence models, 
Moo = 0.90, a = 3.97, Re c = 17.6 x 10 6 . 
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Fig. 8 Comparison of unsteady pressures with experiment, — 0.88, a m = 3.04, d = 0.48, Re c = 17.2 x 10 6 , 
k= 0.6059. 
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Fig. 11 Comparison of unsteady pressures with experiment, M <» = 0.90, a m = 3.97, a = 0.46, Re c = 17.6 x 10 6 , 
k= 0.5919. 
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Fig. 12 Comparison of unsteady pressures between the standard and modified Baldwin- Lomax turbulence models, 
Moo = 0.90, a m = 3.97, & - 0.46, Re c = 17.6 x 10®, k= 0.5919. 
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Clipped Delta Wing 



Fig. 14 Comparison of sectional lift responses between rigid and flexible wings, = 0.90 4° ramping up 
Re c = 15.0 x 10 6 , A= 0.04. 



Fig. 15 Unsteady upper surface pressure responses of rigid wing, M = 0.90, 4° ramping up, Re c = 15.0 x 10 6 , 
A= 0.04. a) 34% section; b) 54% section; c) 68% section; d) 90% section. 
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Fig. 18 Deformation of flexible wing in comparison to rigid wing at 800 time steps, = 0.90, 4° ramping up 
Re c = 15.0 x 10 6 , A= 0.04. 6 



Fig. 19 Responses of sectional lift and elastic angle of attack for flexible wing, M* = 0.90. 4° ramoine ud 
Re c = 15.0 x 10 6 , A= 0.04. 6 
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Clipped Delta Wing (Flexible) 
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Fig. 26 Instantaneous negative- ti- velocity contour plots 
at 68% semispan section, Moo — 0.90, 10° ramping up, 
Re c = 15.0 x 10 6 , A= 0.04. 


Fig. 27 Streamline pattern over the upper surface of 
flexible wing, Moo — 0.90, 10° ramping up, Re c = 
15.0 x 10 6 , A= 0.04. a) 1500 time steps; b) 2400 time 
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Abstract of a proposed paper for the presentation at AIAA Atmospheric Flight Me- 
chanics Conference, August 19-21, 1991, New Orleans, Louisiana. 

Unsteady Navier-Stokes Computations on a Wing-Body 
Configuration in Ramp Motions 

Shigeru Obayashi,* Guru P. Guruswamy** and Eugene L. Tu*** 

NASA Ames Research Center, Moffett Field, California 

Summary 

Unsteady Navier-Stokes computations have been carried out for simulating tran- 
sonic flows over a wing-body configuration undergoing prescribed pitching and ramp 
motions, including flexible modes. The unsteady Navier-Stokes code used in this study 
has been validated by comparing the numerical solutions with experimental data for 
pitching motion cases about wings in transonic flow-fields from a small angle-of-attack 
to a relatively high angle-of-attack, where a leading-edge vortex develops. In this work, 
the method has been extended to handle wing-body configurations. Steady and un- 
steady results at transonic Mach numbers will be shown for the wing-body model and 
comparisons are made with the available experimental and numerical data. The ramp 
motion covers angles of attack where a vortex breakdown is observed in the experiment. 
The effects of unsteadiness will be discussed in detail. 

Introduction 

In the last two decades, there have been extensive developments in computa- 
tional aerodynamics, which constitutes a major part of the general area of computa- 

* Research Scientist, MCAT Institute, San Jose, California. Member AIAA. 

** Research Scientist. AIAA Associate Fellow. 

*** Research Scientist. Member AIAA. 



viscosity is carefully evaluated to capture the leading-edge separation. 

In this abstract, the steady and unsteady results for the ramp motion are pre- 
sented for the rigid wing-body model. In addition, code validation results from previous 
studies 2,4,7 are shown here. Further detailed results including forced flexibility will be 
presented in the full paper. All of the results to be presented in the full paper are new 
and have not been submitted to any other meetings or publications. 

Governing Aerodynamic Equations and Approximations 

The governing equations used in this study is the thin-layer Navier-Stokes equa- 
tions written in conservation-law form in a generalized body-conforming curvilinear 
coordinate system. Several numerical schemes have been developed to solve these equa- 
tions. The present code has two different schemes; the central-difference and upwind 
schemes. Detailed formulas will be given in the full paper. 

The viscosity coefficient in the thin-layer viscous term is computed as the sum of 
Hi + pit where the laminar viscosity, /xj, is taken from the freestream laminar viscosity, 
assumed to be constant for transonic flows. The turbulent viscosity, pi t , is originally 
evaluated by the Baldwin-Lomax algebraic eddy-viscosity model . 10 The modification 
of the turbulence model originally developed for crossflow separation by Degani and 
Schiff 11 is applied to evaluate the turbulent viscosity of flow- fields containing a leading- 
edge vortex. However, for this geometry, this modification has difficulty in locating the 
boundary-layer edge, especially at low angles-of-attack where the leading-edge vortex 
lies near the wing surface. Thus, the Johnson-King model 12 ’ 13 is also applied for 
comparison and the results will be discussed in the full paper. 

Results 


1. Code Calibration on Wings 
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pressures at the outboard section show that the peaks near the leading edge in the 
computed profiles are located more downstream than in the experimental data. This 
indicates that the computed leading-edge vortex location is more inboard than was 
found in the experiment. This discrepancy between computed and measured data is 
widely seen in other numerical results including Ref. 15. Although the Degani-Schiff 
correction 11 is used here, there are problems such as the lack of transition modeling in 
the Baldwin-Lomax model. Nevertheless, the computed results show good qualitative 
agreement with the experiment even for the unsteady cases. 

Next, the ramp motion of the same clipped delta wing was investigated. 4 Figure 2 
shows pressure histories at Moo = 0.9 for the wing ramping up to 5 degrees with a pitch 
rate of A = 0.01. The pitch rate, A, is defined as ac/Uoo. The unsteady computations 
are started from the converged steady-state solution at a = 0.0°. The results show the 
dynamics of vortex and shock development. The corresponding sectional lift curves are 
plotted in Fig. 3. Stall can be seen at 90% semispan section. 

The third example is a comparison of ramp motions of a typical fighter wing 
with different pitch rates. 2 It is well known that the dynamic lift for airfoils can be 
increased by increasing pitch rates. Because dynamic lift is an unsteady phenomenon, 
and is associated with the presence of vortices, it is important to model it accurately. 
In this computation, the Mach number is subsonic and the angle-of-attack varies from 
0° to 30°. Figure 4 shows contour plots of the velocity magnitude at various spanwise 
locations. The plots indicate the presence of the leading-edge vortex. 

To see the effects of increasing pitch rates, computations were made for the three 
pitch rates of 0.025, 0.05 and 0.1. The unsteady sectional lift coefficients are plotted 
against time in Fig. 5. For all three pitch rates, the stall occurs at the wing-tip before 
it occurs at the root section. These plots confirm the increase in the dynamic lift at 
higher pitch rates for this fighter wing case except for the inboard sections. 

Results in this section verify the validation of the present code to accurately 


5 



of-attack is shown in Fig. 7. Comparisons with the experimental data are given at 
wing-span stations of 35, 45, 55 and 85%. The results generally show good agreement 
except at the 85% wing-span where the discrepancies are most likely due to the modeling 
of the wing-tip. Since detailed geometry data was unavailable for the wing-tip of the 
wind-tunnel model, the tip is assumed to be rounded in this study. 

At the inboard spanwise locations, the suction peaks due to the leading-edge 
vortex appear to be slightly underpredicted. This discrepancy becomes worse as angle- 
of-attack increases. However, grid refinement was found to improve the prediction 
significantly. 7 In addition, the turbulence model affects the formation of the leading- 
edge vortex. The performance of the Johnson-King model 12 ’ 13 is currently being inves- 
tigated. In the full paper, these additional numerical issues will be addressed. 

c. Unsteady Results 

In this abstract, the ramp motion case is shown as a sample result. The compu- 
tation was started from the steady-state solution at 4.09° angle-of-attack. Then, the 
model was ramped up 5 degrees with a pitch rate of A = 0.1. After the angle-of-attack 
reached 9.09°, the model was held stationary. Figure 8 shows the grid and pressure con- 
tours on the upper surface of the wing-body model at various angles-of-attack. Figures 
8. a thru 8.d show the grid, pressure contours at a = 4.09° (steady-state), a = 7.59° 
(unsteady) and a = 9.09° (unsteady at maximum a), respectively. 

Figure 9 shows a time history of pressure responses at the 35, 55 and 85% wing 
span locations. The plots show the dynamics of the leading-edge vortex development. 
The corresponding total lift curve during ramp-up is plotted in Fig. 10. The dynamic 
lift increase is clearly shown. After reaching the maximum angle-of-attack, the lift 
coefficient is observed to decrease slowly to the steady-state value. A complete ramp 
motion from 0 to 12 degrees is being performed and will be presented in the full paper. 
In the steady-state results, 7,16 a vortex breakdown was observed at about 12° angle- 
of-attack. The effects of unsteadiness will be discussed in detail. 
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Fig. 6 Close-up view of wing-body grid. 
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b) a = 4.09°, steady-state solution. 



c) a = 7.59°, unsteady solution with pitch rate 0.1. 



d) a = 9.09°, unsteady solution at maximum angle-of-attack with pitch rate 0.1. 


Fig. 8 Instantaneous pressure contours on the upper surface of a wing-body model in 
ramp motion from 4° to 9° angle-of-attack. 
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